Compact Objects for Everyone I: White Dwarf Stars 



CB. Jackson, J. Taruna, S.L. Pouliot, B.W. Ellison, D.D. Lee, and J. Piekarewicz 1 

1 Department of Physics, Florida State University, Tallahassee, FL 32306, USA 

(Dated: February 2, 2008) 

Abstract 

Based upon previous discussions on the structure of compact stars geared towards undergraduate 
physics students, a real experiment involving two upperdevel undergraduate physics students, a 
beginning physics graduate, and two advanced graduate students was conducted. A recent addition 
to the physics curriculum at Florida State University, The Physics of Stars, sparked quite a few 
students' interests in the subject matter involving stellar structure. This, coupled with Stars 
and Statistical Physics by Balian and Blaizot 1] and Neutron Stars for Undergraduates by Silbar 
and Reddy 0], is the cornerstone of this small research group who tackled solving the structure 
equations for compact objects in the Summer of 2004. Through the use of a simple finite-difference 
algorithm coupled to Microsoft Excel and Maple, solutions to the equations for stellar structure 
are presented in the Newtonian regime appropriate to the physics of white dwarf stars. 

PACS numbers: 
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I. OVERVIEW OF THE PROJECT 



It is the central tenet of the present "experiment" that advances in both algorithms and 
computer architecture bring the once-challenging problem of the structure of compact ob- 
jects within the reach of beginning undergraduate students — and even high-school students. 
After a brief historical review in SeclHl a synopsis of stellar evolution is presented in Sec. 11111 
that culminates with a detailed description of the physics of collapsed stars. Following this 
background information, the equations of hydrostatic equilibrium for Newtonian stars are 
derived in Sec. IIVI In particular, the need for an equation of state and the enormous advan- 
tage of scaling the equations are emphasized. The section concludes with the presentation 
of mass-vs-radius relationships for white dwarf stars obtained using both Excel and Maple. 
The special role played by special relativity for the existence of a limiting mass, the Chan- 
drasekhar limit, is strongly emphasized. Finally, a summary and some concluding remarks 
are offered in Sec. I VII 

II. HISTORICAL PERSPECTIVE 

The story of collapsed stars starts in earnest with Subrahmanyan Chandrasekhar in the 
early 1930 's. As a young man of 20, he was embarking from his native India to Cambridge 
University to start life as a graduate student under Fowler's supervision. By 1926 Fowler 
had already explained the structure of white dwarf stars by using the electron degenerate 
pressure only a few months after the formulation of the Fermi-Dirac statistics [3j]. How- 
ever, Chandrasekhar noticed a critical ingredient missing from Fowler's analysis: special 
relativity |4j. Chandrasekhar discovered that as the density of the star increases and the 
momentum of the electrons p becomes comparable with and later exceeds mc, where m is 
the rest mass and c the speed of light, then their ability to support the star against gravity 
weakens. Chandrasekhar concluded that stars with masses above Mch = 1.44M (M a 
solar mass) cannot cool down but will continue to contract and heat up. This limiting mass 
Mch is fittingly known as the Chandrasekhar mass. However, it was not all smooth sail- 
ing for Chandrasekhar. One pre-eminent figure in the field, Sir Arthur Eddington, opposed 
Chandrasekhar publicly and privately with great vigor, even claiming — by arguments found 
generally difficult to follow — that the non-relativistic pressure- density relation should be 
used at all densities. 

Why wouldn't Chandrasekhar silence his critics by revealing the ultimate fate of a heavy 
(M > Mch) star? After all, we now know that such a star will collapse into a neutron 
star, or if too heavy into a black hole. Unfortunately for Chandrasekhar, at the time of his 
ground-breaking discovery it was impossible for him (or for anyone else for that matter) to 
have predicted the existence of neutron stars. It was one year later in 1932 that Chadwick 
proved the existence of the neutron [Hf. From that point on things developed very quickly, 
culminating with the 1933 proposal by Baade and Zwicky that supernovae are created by the 
collapse of a "normal" star to form a neutron star Chandrasekhar was thus vindicated 
and awarded the Nobel Prize in Physics in 1983 for his lifetime contributions to the physical 
processes of importance to the structure and evolution of the stars. 

The process of discovery of these two classes of compact objects (white dwarf and neutron 
stars) is radically different. In the case of white dwarf stars, observation predated Fowler's 
theoretical explanation by more than 10 years. Indeed, in 1915 at the Mount Wilson Ob- 
servatory in California a group of astronomers headed by Walter Sydney Adams discovered 
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that Sirius B — the companion of the brightest star in the night sky — was a white-dwarf 
star, the first to be discovered. On the other hand, it took more than 30 years since the 
bold prediction by Baade and Zwicky p to discover neutron stars. As is often the case in 
science, this discovery marks one of the greatest examples of serendipity. Antony Hewish 
and his team at the Cavendish Laboratory built a radio telescope to study some of the 
most energetic quasi-stellar objects in the universe (quasars). Among the members of the 
team was a young research student by the name of Jocelyn Bell [7J. Shortly after the tele- 
scope started gathering data in 1967, Bell observed a signal of seemingly unknown origin (a 
"scruff"). Particularly puzzling was that the signal showed at a remarkably precise pulse 
rate of 1.33 seconds 8]. So unexpected was the signal, that after a month of futile attempts 
at understanding it, it was dubbed the "Little Green Men." However, by the beginning of 
1968 Hewish, Bell, and collaborators had found three additional pulsating sources of radio 
waves, or pulsars 0|. The final explanation of these enigmatic sources was due to Gold 
shortly after Hewish and Bell published their findings |l0j, [TjJ. Gold suggested that the 
radio signals were due to rapidly rotating neutron stars that rather than emitting pulses of 
radiation, could emit a steady radio signal that it swept around in circles. When the pulsar 
"lighthouse" was pointing in the direction of the telescope, the signal will indeed show up 
as the short "pulses" that Bell had discovered. It is not our intention to present here a 
comprehensive review of either the history or the fascinating phenomena behind pulsars. 
For a recent account see Ref. Li 



III. STELLAR EVOLUTION A SYNOPSIS 

Stellar objects are dynamical systems involving a symbiotic relationship between matter 
and radiation creating enough pressure to oppose gravitational contraction. Thermonuclear 
fusion, the thermally-induced combining of nuclei as they tunnel through the Coulomb 
barrier, is initially responsible for supporting stars against gravitational contraction. The 
ultimate fate of the star depends upon its remaining mass once thermonuclear fusion can 
no longer provide the pressure required to counteract gravity. 

Thermonuclear fusion drives stars through many stages of combustion; the hot center of 
the star allows hydrogen to fuse into helium. Once the core has burned all available hydrogen, 
it will contract until another source of support becomes available. As the core contracts and 
heats, transforming gravitational energy into kinetic (or thermal) energy, the burning of the 
helium ashes begins. For stars to burn heavier elements, higher temperatures are necessary 
to overcome the increasing Coulomb repulsion and allow fusion through quantum-mechanical 
tunneling. Thermonuclear burning continues until the formation of an iron core. Once iron 
- the most stable of nuclei — is reached, fusion becomes an endothermic process. However, 
combustion to iron is only possible for the most massive of stars. When thermonuclear 
fusion can no longer support the star against gravitational collapse, either because they are 
not massive enough (like our Sun) or because they have developed an iron core, the star 
dies and a compact object is ultimately formed. The three final possible stages a star can 
take is a white dwarf, a neutron star, or a black hole. In this our first contribution - 
one that should be accessible to motivated high-school students — we focus exclusively on 
the physics of white dwarf stars. The fascinating topic of neutron stars requires the use of 
general relativity and is therefore reserved for a more advanced forthcoming publication. 

Our Sun will die as a white dwarf star once all of the hydrogen and helium in the core 
has been burned. Towards the final stages of burning, the star will expand and expel most 
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of the outer matter to create a planetary nebula. At the beginning, the non-degenerate 
core contracts and heats up through conversion of gravitational energy into thermal kinetic 
energy. However, at some point the Fermi pressure of the degenerate electrons begins to 
dominate, the contraction is slowed up, and the core becomes a compact object known as 
a white dwarf, cooling steadily towards the ultimate cold, dark, static black dwarf state. 
On the other hand, neutron stars result from one of the most cataclysmic events in the 
universe, the death of a star with a mass much greater than that of our Sun. Electrons in 
these stars behave ultra-relativistically, and as pointed out by Chandrasekhar, hydrostatic 
equilibrium as a cold body becomes impossible to achieve when M>Mch- However, during 
the collapse of the core, a supernovae shock develops ejecting most of the mass of the star 
into the interstellar space and leaving behind an extremely dense core — the neutron star. 
As the star collapses, it becomes energetically favorable for electrons to be captured by 
protons, making neutrons and neutrinos. The neutrinos carry away 99% of the gravitational 
binding energy of the compact object, leaving neutrons behind to support the star against 
further collapse. The pressure provided by the degenerate neutrons, like degenerate electron 
pressure for white dwarf stars, has a limit on the mass it can bear. Beyond this limiting 
mass, no source of pressure exists that can prevent gravitational contraction. If such is the 
case, then the star will continue to collapse into an object of zero radius: a black hole. There 
is a large number of excellent textbooks on the birth, life, and death of stars. The following 
are some references used in this work 0, 

IH HI EE HI Hi- 



A. The Physics of Collapsed Stars 

It is a remarkable fact that quantum mechanics and special relativity — both theories 
perceived as of the very small — play such a crucial role in the dynamics of stars, the former 
in preventing low-mass stars from collapsing into black holes, while the latter in driving the 
collapse of massive stars. In this research experiment we limit ourselves to the study of 
white dwarf (or Newtonian) stars. Further, we assume that white dwarf stars are "cold", 
spherically symmetric, non-rotating objects in hydrostatic equilibrium. 

In dealing with white dwarf stars, the system may be approximated as a plasma containing 
positively charged nuclei and electrons, with the nuclei providing (almost) all the mass and 
none of the pressure and the electrons providing all the pressure and none of the mass. 
This state of matter corresponds to a gas that is electrically neutral on a global scale, but 
locally composed of positively charged ions (nuclei) and the negatively charged electrons. 
Note that even in a zero-temperature, black dwarf state, the matter is effectively ionized. 
For a free atom, the Heisenberg Uncertainty Relation ensures that in the state of minimum 
total energy — kinetic plus electric potential energy — the electrons occupy a finite volume 
of dimension given by the Fermi-Thomas distance, the analogue for more massive atoms 
of the Bohr radius for hydrogen. The density of a collapsed star is so high that the mean 
distance between nuclei is less than the Fermi-Thomas distance: the matter is 'pressure- 
ionized', with the electrons forming an effectively free gas. However, all identical fermions 
(such as electrons, protons and neutrons) obey Fermi-Dirac statistics in that the occupation 
of states is governed by the Pauli Exclusion Principle — no two fermions can exist in the 
same quantum state. For a zero temperature Fermi gas, all available electron states below 
the Fermi energy are filled, while the rest are empty. The Fermi energy is determined solely 
by the electron number density and rest mass. For a compact object such as a white dwarf 
star, the number density is very high and so is the Fermi energy. Typical electron Fermi 
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energies are of the order of 1 MeV which correspond to a Fermi temperature of Tp~10 10 K. 
As the temperature of the system is increased (say from T = to T = 10 6 K or more, as 
estimated for white dwarf interiors), electrons try to jump to a state higher in energy by an 
amount of the order of k^T but fail, as most of these transitions are Pauli blocked. Only 
those high-energy electrons that are within k^T^ 100 eV from the Fermi surface can make 
the transition, but those represent a tiny fraction (T/Tp~10 -4 ) of the electrons in the star. 
Hence, for the purpose of computing the pressure of the system, it is extremely accurate - 
to 1 part in 10 4 — to describe the electrons as a Fermi gas at zero temperature. 



IV. NEWTONIAN STARS 

Let us start by addressing Newtonian stars. For these stars we assume that corrections 
due to Einstein's greatest triumph — the Theory of General Relativity — may be safely 
ignored. White dwarf stars, with escape velocities of only 3% of the speed of light, fall into 
this category. Not so neutron stars — typical escape velocities of half the speed of light 
cause extreme sensitivity to these corrections. 

We start by considering the radial force acting on a small mass element (Am = p(r)AV) 
located at a distance r from the center of the star (see Fig. Q): 

GM(r)Am A w , , A A d 2 r 

'■ - - P(r + Ar)AA + P(r)AA = Am— . (1) 



r 2 " ' v ' dt 2 

Here p(r) is the mass density of the star, M(r) denotes the enclosed mass within a radius 
r, and P is the pressure. Expanding the above equation to lowest order in Ar one obtains 

GM(r)p(r) dP d 2 r 



r 



2 



p(r)** ■ ( 2 ) 



dr dt 2 



Assuming hydrostatic equilibrium (r = f = 0), one arrives at the fundamental equations 
describing the structure of Newtonian stars. That is, 

^ = -^P, P (r =0),P o ; (3a) 
dr r z 

dM „ . . . . . 

— = +Anr 2 p (r) , M r = = , 3b 
dr 

where Eq. (J3b| defines the enclosed mass. 

It is simple to see that in hydrostatic equilibrium, the pressure of the star is a decreasing 
(or at least not increasing) function of r; otherwise the star collapses. Note that the radius 
of the star R is defined as the value of r at which the pressure goes to zero, i.e., P{R) = 0. 
Similarly, the mass of the star corresponds to the value of the enclosed mass at r — R, when 
M = M(R). 



A. Equation of State 

The above set of equations, together with their associated boundary conditions, must be 
completed by an equation of state (EoS), namely a relation P = P(p) between the density 
and pressure. For simplicity, we limit ourselves to the EoS of a zero-temperature Fermi 
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FIG. 1: The radial force acting on a small mass element a distance r from the center of the star. 

gas composed of constituents (e.g., electrons or neutrons) having a rest mass m. The main 
assumption behind the Fermi gas hypothesis is that no correlations (or interactions) are 
relevant to the system other than those generated by the Pauli exclusion principle. For 
some standard references on the equation of state of a free Fermi gas — at both zero and 
finite temperatures — see Refs. |l9l |20||. 

To start, the Fermi wavenumber kp is defined; kp represents the momentum of the fastest 
moving fermion and is solely determined by the number density (n = N/V) of the system, 
where N is the total number of particles in our system, and V is the enclosed volume. 

That is, 

N = 2^&(k F -\k\) = 2 f J^d*ke(k F -\\L\) = V^- 2 , (4) 

or equivalently 

k F = (37rV) 1/3 . (5) 

In Eq. Pj). Q(x) represents the Heaviside (or step) function. Having defined the Fermi 
wavenumber kp, the energy density of the system is obtained from a configuration in which 
all single-particle momentum states are progressively filled in accordance with the Pauli 
exclusion principle. For a degenerate (spin- 1/2) Fermi gas at zero temperature, exactly 
two fermions occupy each single-particle state below the Fermi momentum pp = hkp; all 
remaining states above the Fermi momentum are empty. In this manner we obtain the 
following expression for the energy density: 

^6(A*-|k|)e(k), (6) 

where e(k) is the single-particle energy of a fermion with momentum k. In what follows, the 
most general free-particle dispersion (energy vs. momentum) relation is assumed, namely, 
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one consistent with the postulates of special relativity. That is, 

/ hkc\ 

e(k) = v/ (Hkc) 2 + (mc 2 ) 2 = m<?\J\ + x 2 , I with x = — - I 

\ mc 2 J 



(?) 



In spite of its slightly intimidating form, the integral in Eq. © may be performed in closed 
form. We obtain 

£ = S S{x F ) , (8) 
where So is a dimensionful constant that may be written using dimensional analysis 

\ 1 1 i t "~ ' 

So 



[mc 



(he) 3 ' 

and £{ dimensionless function of the single variable x F = hk F c/mc 2 given by 



(9) 



£(x F ) = — 



x 2 Vl + x 2 dx 



8tt 2 



(1 + 2x 2 ) Jl + x 2 - In (x F + Jl + x 



(10) 

The pressure of the system may now be directly obtained from the energy density by using 
the following thermodynamic relation — which is only valid at zero temperature: 



P 



dVJ 



N,T=0 



PoP- 



(11) 



N,T=0 



In analogy to the energy density, dimensionful and dimensionless quantities for the pressure 
have been defined: 



Pq — £o 
P(xJ = 



[mc 



2\4 



(hef 



(12a) 
(12b) 



It may be surprising to find that a gas of particles at zero-temperature may still generate a 
non-zero pressure. It is quantum statistics, in the form of the Pauli exclusion principle — not 
temperature — that is responsible for generating the pressure. It is nevertheless surprising 
that quantum pressure, a purely microscopic phenomenon, should be ultimately responsible 
for supporting compact stars against gravitational collapse. 

With an expression for the pressure in hand, we are finally in a position to compute its 
derivative with respect to x F (a quantity that we label as rf). As we shall see in the next 
section, rj — a function closely related to the zero-temperature incompressibility — is the 
only property of the EoS that Newtonian stars are sensitive to We obtain 



dP 



Po 



Pn 



x 



3tt 2 



(13) 



x. 



The above expression has a surprisingly simple form that depends on the energy density 
only through its derivatives. Alternatively, one could have bypassed the above derivation in 
favor of the following general relation valid for a zero-temperature Fermi gas: 



dP 

dx-r, 



n- 



de F 

dX-r 



(14) 
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In view of Eq. (|14p. the attentive reader may be asking why go through the trouble of 
computing the energy density and the corresponding pressure if all that is required is the 
dependence of the Fermi energy on x F . The answer is general relativity. While Newtonian 
stars depend exclusively on rj, the structure of relativistic stars (such as neutron stars) are 
highly sensitive to corrections from general relativity. These corrections depend on both the 
energy density and the pressure and will be treated in detail in a future publication. 



B. Toy Model of White Dwarf Stars 

Before attempting a numerical solution to the equations of hydrostatic equilibrium, we 
consider as a warm-up exercise a toy model of a white dwarf star [21] . Assume a white dwarf 
star with a uniform, spherically symmetric mass distribution of the form 

( Po = 3M/4nR\ if r<R; 
HK ' [0 , if r>fl, K ' 

where M and R are the mass and radius of the star, respectively. For such a spherically 
symmetric star, the gravitational energy released during the process of "building" the star 
is given by 

E G = -AnG[ M(r)p(r)rdr , (16a) 
Jo 

rr 

M(r)= / 4nr' 2 p(r')dr' . (16b) 
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For a uniform density star as assumed in Eq. ([15)1. it is straightforward to perform the above 
two integrals. Thus, the gravitational energy released in "building" such a star is given by 

E G (M,R) = - 3 -^f. (17) 

From Eq. (fI7j). we conclude that without a source of gravitational support, a star with a 
fixed mass M will minimize its energy by collapsing into an object of zero radius, namely, 
into a black hole. We know, however, that white dwarf stars are supported by the quantum- 
mechanical pressure from its degenerate electrons, which (at temperatures of about 10 6 K) 
are fully ionized in the star (recall that 1 eV~10 4 K). In what follows, we assume that elec- 
trons provide all the pressure support of the star but none of its mass, while nuclei (e.g., 4 He, 
12 C, . . . ) provide all the mass but none of the pressure. The electronic contribution to the 
mass of the star is inconsequential, as the ratio of electron to nucleon mass is approximately 
equal to 1 : 2000. 

The energy of a degenerate electron gas was computed in the previous section. Using 
Eqs. © and (JBJ) we obtain, 

E F (M, R) = 37r 2 iVm e c 2 ^^, (18) 

where m e is the rest mass of the electron. Naturally, the above expression depends on the 
mass and the radius of the star, although this dependence is implicit in x F . While the toy- 
problem at hand is instructive of the simple, yet subtle, physics that is displayed in compact 
stars, it also serves as a useful framework to illustrate how to scale the equations. 



1. Scaling the Equations 



One of the great challenges in astrophysics, and the physics of compact stars is certainly 
no exception, is the enormous range of scales that one must simultaneously address. For 
example, in the case of white dwarf stars it is the pressure generated by the degenerate 
electrons (constituents with a mass of m e = 9.110 x 10~ 31 kg) that must support stars with 
masses comparable to that of the Sun (M = 1.989 xlO 30 kg). This represents a disparity in 
masses of 60 orders of magnitude! Without properly scaling the equations, there is no hope 
of dealing with this problem with a computer. 

We start by defining fF = Ep/Nm e c 2 from Eq. (fTHj). a quantity that is both dimensionless 
and intensive (i.e., independent of the size of the system). That is, 

UM = 3- 2 ^, (19) 

where a closed-form expression for S(x F ) has been displayed in Eq. (fTT)|) . Note that the 
scaled Fermi momentum x F quantifies the importance of relativistic effects. At low density 
(x F <C 1) the corrections from special relativity are negligible and electrons behave as a 
non- relativistic Fermi gas. In the opposite high-density limit (x F ^>l) the system becomes 
ultra-relativistic and the "small" (relative to the Fermi momentum) electron mass may be 
neglected. We shall see that in the case of white dwarf stars, the most interesting physics 
occurs in the x F ~ 1 regime. 

The dynamics of the star consists of a tug-of-war between gravity that favors the collapse 
of the star and electron-degeneracy pressure that opposes the collapse. To efficiently com- 
pare these two contributions, the contribution from gravity to the energy must be scaled 
accordingly. Thus, in analogy to Eq. (|19|) . we form the corresponding dimensionless and 
intensive quantity for the gravitational energy (f G = E G /Nm e c 2 ): 

, ^ 3 ( GM\ ( M \ 3 (GM\ ( m n \ 

where we have assumed that the mass of the star, M = Am n , may be written exclusively in 
terms of its baryon number A and the nucleon mass m n (the small difference between proton 
and neutron masses is neglected). This is an accurate approximation as both nuclear and 
gravitational binding energies per nucleon are small relative to the nucleon mass. Further, 
Y e = Z/A represents the electron-per-baryon fraction of the star (e.g., Y e = l/2 for 4 He and 
12 C, and Y e = 26/56 in the case of 56 Fe). 

The final step in the scaling procedure is to introduce dimensionful mass M and radius 
Rq, quantities that, when chosen wisely, will embody the natural mass and length scales in 
the problem. To this effect we define 

M = M/M and R = R/R . (21) 

In terms of these natural mass and length scales, the gravitational contribution to the energy 
of the system takes the following form: 



While the dependence of the above equation on M and R is already explicit, the Fermi gas 
contribution to the energy depends implicitly on them through x F . To expose explicitly the 
dependence of f F on M and R we perform the following manipulation aided by relations 
derived in Sec. IIV Al 



/ Hk F c\' 
\m e c 2 J 



9tt 



Y 



M \ fhc/m e c 2 
m r 



R 







M 



(23) 



We have already referred earlier to M and Rq as the natural mass and length scales in 
the problem, but their values have yet to be determined. Thus, they are still at our disposal. 
Their values will be fixed by adopting the following choice: let the "complicated" expressions 
enclosed between square brackets in Eqs. (|2*2*|l and (|2*3j) be set equal to one. That is, 



GM C 
R^ 



m r 



Y P m ( 



9^ 



Y, 



M \ fhc/m e c 2 



R 







(24) 



This choice implies the following values for white dwarf stars with an electron-to-baryon 
ratio equal to Y e = 1/2: 



M 

Rq 



-VT&rctG 1 m n Y? = 10.599 M Y e 2 — ► = 2.650 M , 

Ye 1 /2 



2 G 



7T -1/2 

a, 1 



he 



m P cr 



Y P = 17250 kmY P 



' Y e =l/2 



8 623 km . 



(25a) 
(25b) 



Here the minute dimensionless strength of the gravitational coupling between two nucleons 
has been introduced as 



he 



5.906 x 10 



-39 



(26) 



The aim of this toy-model exercise is to find the minimum value of the total (gravitational 
plus Fermi gas) energy of the star as a function of its radius for a fixed value of its mass. 
Before doing so, however, a few comments are in order. First, from merely scaling the 
equations and with no recourse to any dynamical calculation we have established that white 
dwarf stars have masses comparable to that of our Sun but typical radii of only 10 000 km 
(recall that the radius of the Sun is R Q ~ 700 000 km). Further, we observe that while Rq 
scales with the inverse electron mass, the mass scale M is independent of it. This suggests 
that neutron stars, where the neutrons provide all the pressure and all the mass, will also 
have masses comparable to that of the Sun but typical radii of only about 10 km. 

Now that the necessary "scaling" machinery has been developed, we return to our original 
toy-model problem. Taking advantage of the scaling relations, the energy per electron in 
units of the electron rest energy is given by: 



f(M,x F ) = f G (M,x F ) + f F (x F ) = -M 2/3 x F + 37r 2 ^ , r ' ' 



X, 



(27) 



The mass-radius relation of the star may now be obtained by demanding hydrostatic equi 

librium: 

/ df(M,x F ) \ =Q 

V dx F )m 



(28) 
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While a closed-form expression has already been derived for the energy density S(x F ) in 
Eq. (JlUj) . it is instructive to display explicit non-relativistic and ultra- relativistic limits. 
These are given by 



To conclude this section an output from a Maple code has been included to illustrate how 
simple, within the present approximation, it is to compute the radius of an arbitrary mass 
star (see Fig. EJ). For the present example, a 1 M star has been used. First, scales for 
input quantities, such as the dimensionful mass (M ) and length (i?o), are defined. Next, 
energies and their derivatives are computed and a plot displaying the latter is generated. 
The derivative of the gravitational energy (actually the negative of it) is constant and is 
displayed with a black horizontal line. Similarly, the derivative of the Fermi energy in the 
ultra-relativistic limit (horizontal green line) is also a constant equal to 3/4, independent 
of the mass of the star. The blue line with a constant slope displays the derivative of the 
Fermi gas energy in the non-relativistic limit. Finally, the exact Fermi gas expression, which 
interpolates between the non-relativistic and the relativistic result, is displayed with the red 
line. 

The equilibrium density of the star is obtained from the intersection of the red and blue 
lines with the gravitational line. In the non-relativistic case the solution may be computed 

2/3 

analytically to be xpo = 5M /3. However, this non-relativistic prediction overestimates the 
Fermi pressure and consequently also the radius of the star. The non-relativistic predictions 
for the radius of a 1 solar-mass white dwarf star is i?NR = 7162 km. In contrast, the result 
with the correct relativistic dispersion relation is considerably smaller at an .Rnr = 4 968 km. 
Yet an even more dramatic discrepancy emerges among the two models. While the non- 
relativistic result guarantees the existence of an equilibrium radius for any value of the star's 

mass (-Rnr = 3/5M 1//3 ), the correct dispersion relation predicts the existence of an upper 
limit beyond which the pressure from the degenerate electrons can no longer support the star 
against gravitational collapse. This upper mass limit, known as the Chandrasekhar mass, is 
predicted in the simple toy model to be equal to: 



As it will be shown later, accurate numerical results yield (for Y e = 1/2) a Chandrasekhar 
mass of M Ch = 1.44 M Q . Thus, not only does the toy model predict the existence of a 
maximum mass star, but it does so with an 80% accuracy. 

C. Numerical Analysis 

We now return to the exact (numerical) treatment of white dwarf stars. While the toy 
model problem developed earlier provides a particularly simple framework to understand 
the interplay between gravity, quantum mechanics, and special relativity, a quantitative 
description of the systems demands the numerical solution of the hydrostatic equations 
[Eq. (J3J)]. For the present treatment, however, we continue to assume that the equation of 
state is that of a simple degenerate Fermi gas. In this case the equation of state is known 
analytically and it is convenient to incorporate it directly into the differential equation. In 




x F < 1 ; 
x F > 1 . 



(29) 



M Ch = (3/4) 3/2 M = 1.72 M ( 



(30) 
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[ Toy Model of Compact Stars: 



[ > restart : 

> MO := 2 . 651 : 

> RO := 8625: 
[ > M := 1 . OOO : 
[ > Mbar := M/MO : 
[ Compute energies: 

[ > f G := -Mbar^ (2/3) *xF : 

[ > f F := (3/xF^3) *int (x A 2*sqrt (l+x A 2) 

[> f FNR := 1+ (3/10) *xF A 2 : 

[ > f FUR := (3/4) *xF: 

[ Compute derivative of the energies: 

[ > fGl := dif f (fG, xF) : 

[> fFl := simplify (dif f (fF, xF) ) : 

simplify (dif f (fFNR, xF) ) : 
simplify (dif f ( f FUR, xF) ) : 



#dimensionf ull mass scale (in solar masses) 
#dimensionf ull length scale (in km) 
#compute structure for a 1 solar-mass star 
#dimensionless mass 



#gravitat ional energy 
x=0..xF): #Fermi energy 
#Fermi energy (NR limit) 
#Fermi energy (UR limit) 



[ > fFNRl 
[ > fFURl 
[ Make a plot: 
> plot ( [-fGl, fFl, fFNRl, fFURl] , xF = 



#gravitat ional energy 
#Fermi energy 
#Fermi energy (NR limit) 
#Fermi energy (UR limit) 



2, 



co lor= [black, red, blue, green] , thickness=3) ; 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
xF 



[ Compute the radius of a 1 solar mass star: 

[> xFO := f solve (fGl+fFl, xF=l) 
[> RbarO := Mbar A ( 1 / 3 ) /xFO : 
[> Mass := Mbar*M0: 
[> Radius := Rbar0*R0: 



tsolve for the density 
tsolve for the Radius 
textract dimensionf ull mass 
textract dimensionf ull radius 



> printf ("M=%4 .2f MSolar; R=%4.0f km" , Mass , Radius ) ; 
M=1.00 MSolar; R=4968 km 



FIG. 2: Maple code displaying the interplay between various physical effects on the "toy- model' 
problem of a M = 1M & star. in 



this way Eq. (|5a|) becomes 



dx F 
dr 



GM{r)p{r) 



r 2 r) 



(31) 



where the equation of state enters only through a quantity directly related to the zero- 
temperature incompressibility 20]. This quantity, rj = dP/dx F , was defined and evaluated 
in Eq. (JTHJ). Moreover, the density of the system p{r) is easily expressed in terms of x F . It 
is given by 

o2~ 



P 



he 



m r , 



3ir 2 Y P F ' 



(32) 



At this point all necessary relations have been derived and the equations of hydrostatic 
equilibrium, first displayed in Eq. but with an equation of state still missing, may now 
be written in the following form: 



dx F 
dr 

d~M 

dr 



(GMo 
Woe 2 



Y e m f 



M "t/ 1 + X r 



x F {r 



+ 



3^ 



Y 



M \ ( hc/m e c 2x " 



m r 



Rq 



—2 3 

r x F 



= 0)=x Fc ; (33a) 
M(f=0) = 0. (33b) 



Here the dimensionless distance r and the central (scaled) Fermi momentum x Fc have been 
introduced. The structure of the above set of differential equations indicates that our goal 
of turning Eq. (JBJ) into a well-posed problem, by directly incorporating the equation of state 
into the differential equations, has been accomplished. But we have done better. By defining 
the natural mass and length scales of the system (M and Rq) according to Eq. pijl. the two 
long expressions in brackets in the above equations reduce to the simple numerical values of 
5/3 and 1/3, respectively. Finally, then, the equations of hydrostatic equilibrium describing 
the structure of white dwarf stars are given by the following expressions: 



dx F 
dr 

d~M 
dr 



f{r;x F ,M) , x F (r = 0)=x Fc ; 
g(f;x F ,M), M(r = 0)=0, 



(34a) 
(34b) 



where the two functions on the right-hand side of the equations (/ and g) are given by 



f(r; x F ,M) 



5M (/1 



x. 



and g(r;x F ,M) = +3r 2 x 3 



3 r 



tt2 



Xr 



(35) 



These coupled set of first-order differential equations may now be solved using standard nu- 

However, for those students not 
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merical techniques, such as the Runge-Kutta algorithm 
yet comfortable with writing their own source codes, the use of an "off-the-shelf" spreadsheet 
(here Microsoft Excel has been used), together with a crude low-order approximation for the 
derivatives has been shown to be adequate. As in the toy-model problem, solutions will be 
presented using the full relativistic dispersion relation as well as the non-relativistic approx- 
imation, where in the latter case the square- root term appearing in the function /(f; x F , M) 
is set to one. 
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V. NUMERICAL TECHNIQUES FOR EVERYONE 



In this section we present numerical solutions for the structure of Newtonian (white- 
dwarf) stars by employing a variety of numerical techniques and programming tools. 

1 . White Dwarf Stars with Excel 

This section should be ideal for those students with a basic knowledge of calculus and with 
no programming skills. High school students that have learned the concept of derivatives in 
their introductory calculus class should be able to complete this part of the project with no 
problem. Indeed, one could use the definition of the derivative of a function F(x) 

x i F(x + h) - Fix) . , 

F'(x) = im -A 1 U. , (36) 

h— >0 II 

to approximate the value of the function at a neighboring point x + h. That is, 

F{x + h) = F{x) + hF\x) + 0(h 2 ) . (37) 

The term 0(h 2 ) indicates that the error that one makes in computing the value of the 
function at the neighboring point scales as the square of h. Thus, for this low-order approx- 
imation of the derivative, we selected a very small value of h in order to ensure numerical 
accuracy. Other higher-order algorithms, such as the venerated Runge-Kutta method, con- 
tain errors that scale as (9(/i 4 ). These algorithms can attain the same degree of numerical 
accuracy as the one presented here with a dramatic reduction in computational time. 

So how do we turn Eq. ()37|) to our advantage? By simply looking at the structure of the 
(scaled) equations of hydrostatic equilibrium Eq. (|34j) we can readily write 

x F (f + Ar) = x F (f) + Ar f(f; x F , M) , (38a) 
M(r + Ar) = M(r) + Ar g(r; x F ,M). (38b) 

The resulting difference equations are recursion relations that enables one to "leapfrog" 
from point to point in a grid for which the various points are separated by a fixed distance 
Af. Recursion relations such as this one are particularly well suited to be solved with a 
spreadsheet. Fig. El shows the results produced using Excel (lines). 

To start the solution of the difference equations one notes that the right-hand side of the 
two recursion relations given above are completely known at r = (recall that appropriate 
boundary conditions have already been specified). This enabled us to compute both the 
scaled Fermi momentum and the enclosed mass on the next grid point r = Ar. With this 
knowledge we could again evaluate the right-hand side of the recursion relations but now at 
r = Ar. Now the values of x F and M at the next grid point (r = 2Ar) may be computed. We 
continued in this manner until the Fermi momentum goes to zero. This point defines the 
radius of the star, while the mass of the star is the value of the enclosed mass at this last 
point (this point also corresponds to when the pressure goes to zero). Up to this point, both 
radius and mass are obtained in dimensionless units. To convert back to physical units we 
simply multiplied these dimensionless quantities by the dimensionful parameters (Rq and 
Mq) defined in Eq. (J25j) . Of course, there is no need to repeat the full calculation for a 
different value of Y e as we scaled Rq and M appropriately. In order to create the complete 
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mass-radius relation, we repeated the same procedure for a large number of central densities. 
Here a step size of Ar = 0.0001 was used throughout and scaled central densities ranged from 
0.1 to 100, moving in steps of 0.1 at first, then to increments of 1. Doing so produces Fig. El 
and, in particular, a Chandrasekhar limit very close to 1.4 M Q . However, as alluded in the 
toy-model problem and confirmed in this numerical calculation, there is no Chandrasekhar 
limit if one uses a non-relativistic dispersion relation. 




FIG. 3: Mass- us-radius relation for white dwarf stars obtained using Excel (lines) and Maple 
(symbols). 



2. White Dwarf Stars with Maple 

The Maple code was created in a simple manner. It starts with a switch (or option) that 
queries the user for the type of Fermi gas equation of state (relativistic or non-relativistic) 
to be used. With this information one uses Eq. (J23J) to define appropriate dimensionful mass 
and radius parameters which are used after the scaled equations have been solved. Next, 
a Maple procedure was written to solve the scaled differential equations (consistent with 
the switch) for a given central (scaled) Fermi momentum using a classical numeric method. 
We noted that since in the structure equations the radius appears in the denominator of 
several expressions, we could not start calculating at a zero radius. Of course, as the limit 
is well defined and finite as the radius goes to zero, one could use an extremely small value 
to solve this problem (we used a scaled radius of 10~ 23 for the first point). For a given 
central Fermi momentum, the equations were numerically solved for x F and the enclosed 
mass M as a function of the scaled radius r using a very small step size of Ar = 0.00025. 
Such a small value for Ar is necessary in order to account for the rapidly- varying behavior 
of the density on the surface of the star. The program was instructed to stop once the scaled 
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Fermi momentum turned negative, with the previous point defining the radius of the star 
and the value of the enclosed mass at this point representing the mass of the star. Physical 
dimensions were restored by multiplying these scaled values by the dimensionful mass and 
radius parameters computed earlier. Having done this once, the procedure was repeated for 
a large range of central densities so one could accurately map the mass-i>s-radius relation of 
the star. Once these values were stored, a plot was generated (Fig. EJ). Fig. |H] illustrates the 
consistency of the results using either Excel or Maple; the data calculated using Excel are 
shown as lines, and the data from Maple are shown as symbols. 



VI. SUMMARY - CONCLUDING REMARKS 

Isaac Newton once said: If I have seen farther than others, it is because I was standing on 
the shoulders of giants. The foundations for the present experiment were laid by giants such 
as Chandrasekhar, Fermi, Dirac, and Einstein. Their seminal work placed the fascinating 
world of compact stars within the reach of the whole scientific community. Literally, they 
reduced the problem of compact stellar objects to quadratures. When these insights are 
combined with the remarkable advances in computer processing and algorithms, the prob- 
lem becomes accessible to undergraduate students. (Indeed, highly motivated high-school 
students should also be able to tackle most of this problem.) The exercise reported in this 
manuscript is a testimony to this fact. 

Our work benefited greatly from earlier contributions by Balian and Blaizot [H, and Silbar 
and Reddy that pioneered the idea of bringing the physics of stars to the realm of the 
classroom. Here we have followed closely on their footsteps. For the present project a group 
involving two upper-level undergraduate physics students, a beginning physics graduate, and 
two advanced graduate students was assembled. It was demonstrated that with a limited 
knowledge of calculus and physics, undergraduate students can easily tackle the structure 
equations for white dwarf stars. Moreover, we are convinced that with a relative small 
amount of mentoring, the same will be true for motivated high school students 0. 

Students learned several important lessons from this project. One of them relates to the 
usefulness of scaling the equations-without scaling, the problem would have been unsolvable. 
This is due to the tremendous range of scales encountered in this problem; there are more 
than 60 orders of magnitude between the minute electron mass and the immense solar mass. 
Another important lesson learned is that, contrary to what seems to happen in the classroom, 
most problems in physics have no analytic solution. Thus, numerical analysis is a necessary 
step towards a solution. It was shown that with a limited knowledge of calculus one can 
derive suitable recursion relations to arrive at accurate solutions to the differential equations 
by using a simple spreadsheet program like Microsoft Excel. For more advanced students, 
symbolic programs (such as Maple or Mathematica) provide a more efficient method for 
arriving at the solutions. Due to licensing agreements, Maple was used in this work. 

In conclusion, this project successfully utilized common resources to solve structure equa- 
tions for compact objects in the Newtonian regime. Building on this project, students are 
now in a position to study the fascinating physics of neutron stars. The structure of neutron 
stars, however, poses several additional challenges. First and foremost, Newtonian gravity 
must be replaced by general relativity. This implies that the structure equations must be 
replaced by the Tolman-Oppenheimer-Volkoff equations k| 12, IB- Second, at the higher 



densities encountered in the interior of neutron stars, the equation of state receives impor- 
tant corrections from the interactions among the neutrons. That is, Pauli correlations are 
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no longer sufficient to describe the equation of state and some realistic equations of state 
should be used 0, H2| • This topic of intense research activity is of relevance to the physics 
of neutron stars and to the structure of those exotic compact objects known as hybrid and 
quark stars. 
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